#!/usr/bin/python3
# by Sun Smallwhite <niasw@pku.edu.cn>(https://github.com/niasw)

import numpy
import matplotlib
matplotlib.use('Agg'); # plot and save, never showing

import matplotlib.pyplot as plt

def calcXvector(clusters,dimension,numvec):
  '''
# for a paper p:
#  x(p)=sqrt(|#p|)/sqrt(|#c(p)|)
# input:
#  clusters: the division result
#  numvec: paper number of initial clusters
#  dimension: number of initial clusters
  '''
  xvec=numpy.zeros(dimension)+1.;
  for c in clusters:
   csum=numpy.array(numvec)[c].sum();
   for p in c:
    xvec[p]=numpy.sqrt(numpy.float64(numvec[p])/numpy.float64(csum));
  return xvec;

def monitorXandErr(eigvec,clusters,dimension,numvec):
  '''
# monitoring x vector and relative error
  '''
  xvec=calcXvector(clusters,dimension,numvec);
  reerr=1-xvec.dot(eigvec)**2/xvec.dot(xvec)/eigvec.dot(eigvec);
  return xvec,reerr;

def visualizeCombine(eigvec,clusters,numvec):
  '''
# visualize x and eigvec
  '''
  def _idxfunc(x):
   return x[1];
  def _keyfunc(x):
   return x[0];
  sortmap=sorted([(eigvec[p],p) for p in range(0,len(eigvec))],key=_keyfunc);
  xvec=calcXvector(clusters,len(eigvec),numvec);
  mappedxvec=numpy.array([xvec[_idxfunc(sortmap[p])] for p in range(0,len(eigvec))]);
  mappedeigv=numpy.array([eigvec[_idxfunc(sortmap[p])] for p in range(0,len(eigvec))]);
  print("len eigv: "+str(len(eigvec)));
  print("max eigv: "+str(mappedeigv[len(eigvec)-1]));
  print("min eigv: "+str(mappedeigv[0]));
  print("avg eigv: "+str(mappedeigv.sum()/len(eigvec)));
  print("neg eigv: "+str((mappedeigv<0).sum()));
  maponexvec=mappedxvec/numpy.sqrt(len(clusters));
  plt.clf();
  plt.plot(range(0,len(eigvec)),mappedeigv,color="green",linewidth=2,linestyle='--',label='eigenvector');
  plt.plot(range(0,len(eigvec)),maponexvec,color="blue", linewidth=2,linestyle='-',label='x vector');
  plt.grid();
  plt.xlabel('paper index (sorted by eigvec)');
  plt.ylabel('vector component');
  plt.legend(loc='upper left');
  plt.savefig('../out/xveig'+str(len(clusters))+'.png');
  # plt.show(); # don't forget to commentize "matplotlib.use('Agg');" in the beginning

def clussprs2clusters(clussprs,sortmap):
  '''
# generate clusters set from cluster seperator data
  '''
  def _idxfunc(x):
    return x[1];
  clusters=[list(set([_idxfunc(sortmap[p]) for p in range(clussprs[c],clussprs[c+1])])) for c in range(0,len(clussprs)-1)];
  return clusters;
  
def calcClusters(eigvec, numvec, clusnum=20, monitor=False):
  '''
# combine paper clusters to make x approching eigv
# by Sun Smallwhite <niasw@pku.edu.cn>(https://github.com/niasw)

# input:
#  eigvec: eigenvector of input flow matrix A (indexed by input clusters)
#  numvec: number of papers in each input clusters
#  clusnum: max number of clusters
#  monitor: whether visualize the progress
# output:
#  clusters: clusters set [[paper indices],...]
  '''
  if (len(eigvec)!=len(numvec)):
   raise Exception('Error: dimension not match, eigvec and numvec.');
  dimension=len(eigvec);
  totalpnum=int(numpy.round(numpy.array(numvec).sum()));
  curcnum=dimension;
  setcnum=clusnum;
  bestrcd=numpy.abs(clusnum-curcnum); # best record
  noimprovecount=0; # prevent loop

  def _compare(x,y):
   if (x[1]==y[1]):
    return x[2]-y[2]; #large eigvec should join first
   else:
    return x[1]-y[1];

  def cmp_to_key(mycmp):
   '''
( python 2 to 3 )
Convert a cmp= function into a key= function
http://code.activestate.com/recipes/576653-convert-a-cmp-function-to-a-key-function/
Created by Raymond Hettinger on Wed, 18 Feb 2009 (MIT)
   '''
   class K(object):
    def __init__(self, obj, *args):
     self.obj = obj
    def __lt__(self, other):
     return mycmp(self.obj, other.obj) < 0
    def __gt__(self, other):
     return mycmp(self.obj, other.obj) > 0
    def __eq__(self, other):
     return mycmp(self.obj, other.obj) == 0
    def __le__(self, other):
     return mycmp(self.obj, other.obj) <= 0  
    def __ge__(self, other):
     return mycmp(self.obj, other.obj) >= 0
    def __ne__(self, other):
     return mycmp(self.obj, other.obj) != 0
   return K;

  def _searchless(xs,value):
   for idx in range(0,len(xs)):
    if (xs[idx][1]<=value):
     return idx;
   return -1;

  if (monitor):
   print('1st try:'+'\n  now cluster number: '+str(dimension)+'\n  set cluster number: '+str(setcnum));
  if (dimension<=clusnum):
   print('Setting Cluster Number is smaller than Current Cluster Number, reset to equal.');
   clusnum=dimension;
   setcnum=dimension;
   curcnum=dimension;
   return [[idx] for idx in range(0,len(numvec))];
  while (curcnum!=clusnum and noimprovecount<=7): # no improvement for 7 times => treat as loop
   lastpush=0; # keep the eigvec order undisturbed, last paper number
   lastclus=None; # keep the eigvec order undisturbed, last cluster name
   curcnum=0; # current cluster number
   clusters=[]; # return
   tmpclus=[];
   curpnum=0; # current paper number in current cluster
   midpnum=1; # setting paper number in current cluster
   stacks={};
   for it in range(0,dimension): # group by target number
    stackname=numpy.float64(numvec[it])/eigvec[it]**2/setcnum;
    if (stackname<1.): # too small
     stackname=1;
    elif (stackname>totalpnum): # too large
     stackname=totalpnum;
    else: # normal
     stackname=int(numpy.ceil(stackname));
    if stacks.__contains__(stackname):
     stacks[stackname].append((it,numvec[it],eigvec[it]));
    else:
     stacks[stackname]=[(it,numvec[it],eigvec[it])];
   tmpstacks=[(it,stacks[it]) for it in sorted(list(stacks))]; # sort dict
   stack=[];
   acumu=0; # \int_{0}^{midpnum} targnum(curpnum)\, d curpnum
   while (len(tmpstacks)>0 or len(stack)>0): # there are untouched clusters
    if (midpnum!=curpnum and len(stack)==0): # call next targnum
     stackobj=tmpstacks.pop(0);
     targnum=stackobj[0];
     midpnum=int(numpy.ceil((targnum+numpy.sqrt(targnum**2-(targnum*curpnum-acumu)*4.))*0.5)); # (acumu + targnum*(midpnum-curpnum))/midpnum = midpnum
     stack=sorted(stackobj[1],key=cmp_to_key(_compare),reverse=True); # sort by numvec[it]
    idx=_searchless(stack,midpnum-curpnum); # use the medium targnum (better than latest)
    if (idx==-1): # cluster not found
     if (len(tmpclus)==0): # targnum too small
      for rit in stack: # one cluster complete
       tmpclus.append(rit[0]);
       lastpush=numvec[rit[0]];
       clusters.append(tmpclus);
       tmpclus=[];
       curpnum=0;
       acumu=0;
       curcnum+=1;
      stack=[];
     else: # completed a target cluster
      lastpush=numpy.array(numvec)[tmpclus].sum();
      clusters.append(tmpclus);
      tmpclus=[];
      curpnum=0;
      acumu=0;
      curcnum+=1;
    else: # put the cluster in
     curpnum+=stack[idx][1];
     acumu+=targnum*stack[idx][1];
     tmpclus.append(stack.pop(idx)[0]);
   if (midpnum!=curpnum and curpnum!=0): # lack of clusters
    if (numpy.array(numvec)[tmpclus].sum()>=lastpush): # right order
     clusters.append(tmpclus);
     tmpclus=[];
     curpnum=0;
     acumu=0;
     curcnum+=1; # flush all we have
    else:
     clusters[len(clusters)-1].extend(tmpclus);
   setcnum=int(numpy.ceil(setcnum*numpy.sqrt(clusnum/curcnum))); # setting cluster number
   if (numpy.abs(clusnum-curcnum)<bestrcd):
    bestrcd=numpy.abs(clusnum-curcnum);
    noimprovecount=0;
   else:
    noimprovecount+=1;
   if (monitor):
    print('result :'+'\n  now cluster number: '+str(curcnum)+'\n  set cluster number: '+str(setcnum));
    xvec,reerr=monitorXandErr(eigvec,clusters,dimension,numvec);
    print('unparallel ratio = '+str(reerr));
   if monitor:
    visualizeCombine(eigvec,clusters,numvec);
    
  clusternumvec=numpy.array([sum(numpy.array(numvec)[clusters[it]]) for it in range(0,len(clusters))]);
  sqrtclusternumvec=numpy.sqrt(clusternumvec);  
  eigclussum=numpy.array([sum(numpy.array(eigvec)[clusters[it]]) for it in range(0,len(clusters))]);
  flowdotproduct=sum(eigclussum/sqrtclusternumvec);
  
  clustersizearr=numpy.array([len(clusters[it]) for it in range(0,curcnum)]);
  clustervectorlength=sum(clustersizearr/(sqrtclusternumvec*sqrtclusternumvec));
  return clusters;
